pacman::p_load(tidyverse, data.table, sf, raster)
rm(list=ls())
################################################################################ 


####################### Extract carbon stocks from rasters into geographic unit 

for (shpfile_name in c("argentina_county","brazil_county") ){
  
  if(shpfile_name=="argentina_county") {
    region = 'argentina'
    raster_bb = c(-76,-52, -57,-20)
  } else {
    region = 'brazil'
    raster_bb = c(-76,-34,-35, 6) 
  }
  shpfile = sf::read_sf(dsn=path.expand(c("../../data/shapefiles/clean/")), layer=shpfile_name)
  
  for (carbonstock in c("aboveground","belowground") ){
    
    #Carbon stock raster
    raster_file = crop(raster(paste0("../../data/environmental/raw/",carbonstock,".asc")), extent(raster_bb)) 

    #Extract data
    mean_stock <- data.frame(extract(raster_file,  shpfile, fun = mean, na.rm = TRUE, df=TRUE, sp=TRUE))
    median_stock <- data.frame(extract(raster_file,  shpfile, fun = median, na.rm = TRUE, df=TRUE, sp=TRUE))
    max_stock <- data.frame(extract(raster_file,  shpfile, fun = max, na.rm = TRUE, df=TRUE, sp=TRUE))
    setnames(mean_stock, old = carbonstock , new = 'mean' )
    setnames(median_stock, old = carbonstock , new = 'median' )
    setnames(max_stock, old = carbonstock , new = 'max' )
    
    #Merge statistics into one file
    merge_cols = c('ctry','st','st_id','cty','cty_id')
    carbonstock_df = merge.data.frame(mean_stock, median_stock, by=merge_cols)
    carbonstock_df = merge.data.frame(carbonstock_df, max_stock, by=merge_cols)
    carbonstock_df <- carbonstock_df %>% mutate(across(.cols = everything(), toupper))
    carbonstock_df$type = carbonstock
    
    #Append above and belowground files into one country file
    if(carbonstock == 'aboveground') {
      final_df = carbonstock_df
    } else { 
      final_df = rbind(final_df,carbonstock_df)}
    
  }
  
  write.csv(final_df, gzfile(paste0("../../data/environmental/clean/carbonstocks_", region, "_extracted.csv.gz")), row.names = FALSE)
}


####################### Create combined Argentina+Brazil dataframe

for (region in c('argentina','brazil') ){
  
  df = fread(paste0("../../data/environmental/clean/carbonstocks_",region,"_extracted.csv.gz"))
  df = df[,c('cty_id','mean','median','max','type')]
  df_1 = rename(df, c("county_id"="cty_id"))
  
  df = fread(paste0("../../data/landuse/clean/geographicunits_",region,".csv.gz"))
  df_2 = df[,-c('land_area_km2','population')]
  
  df = merge.data.frame(df_2, df_1, by='county_id')
  write.csv(df, gzfile(paste0("../../data/environmental/clean/carbonstocks_",region,".csv.gz")), row.names = FALSE)
  
}

df_ar = fread("../../data/environmental/clean/carbonstocks_argentina.csv.gz")
df_br = fread("../../data/environmental/clean/carbonstocks_brazil.csv.gz")
df_br = df_br[,-c('microregion','mesoregion','microregion_id','mesoregion_id','amc_id')]
df_arbr = bind(df_ar,df_br)
write.csv(df_arbr, gzfile("../../data/environmental/clean/carbonstocks_argentinabrazil.csv.gz"), row.names = FALSE)


